Interpolate Subroutine

SUBROUTINE Interpolate &
!
(cvobs, cvest, controlpts, est, var, neighbors, sill)
	
IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short),              INTENT(IN) :: neighbors
REAL (KIND = float),                 INTENT(IN) :: sill
REAL (KIND = float), DIMENSION(:),   INTENT(IN) :: cvest
REAL (KIND = float), DIMENSION(:,:), INTENT(IN) :: cvobs
REAL (KIND = float), DIMENSION(:),   INTENT(IN) :: controlpts

!arguments with intent(out):
REAL (KIND = float), INTENT(OUT) :: est !!estimated value
REAL (KIND = float), INTENT(OUT) :: var !!variance of estimated value

!Local variable declarations
INTEGER (KIND = short) :: n
!REAL,DIMENSION(neighbors+1) :: musv
REAL (KIND = float) :: ck
!--------------------------end of declarations---------------------------------

n = SIZE(controlpts)

!Derive Weights
weights = MATMUL(cvobs,cvest)
		
!Check weights sum to 1 (n=1 th element is the lagrange parameter)
ck = SUM(weights(1:neighbors))
IF (NINT(ck) .NE. 1) THEN
    CALL Catch ('warning', 'Kriging interpolation',  &
          'Kriging weights do not sum to 1' )
END IF
			 
!Estimate datum value sum(weights*observations)
est = SUM ( weights(1:neighbors) * controlpts )

!Estimate the local mean: derive new set of weights from the observtions covariance
!matrix and an estimateion covariance vector which is equal to zero. This filters
!the residual component of the estimate, returning the local mean.
!musv=0
!musv(neighbors+1)=1
!musv=MATMUL(cvobs,musv)
!localmu=SUM(musv(1:neighbors)*controlpts)

!Derive Estimation Variance: ssum(weights*semivariances) + lagrangian
var = sill - SUM(weights(1:n)*cvest(1:n)) + weights(n+1)

IF ( var < 0. ) THEN
    CALL Catch ('warning', 'Kriging interpolation',  &
          'Negative variance produced! Check validity &
            of semivariogram model' )
END IF

RETURN
END SUBROUTINE Interpolate